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Abstract 

The present work considers the localization problem in wireless sensor networks 
formed by fixed nodes. Each node seeks to estimate its own position based on 
noisy measurements of the relative distance to other nodes. In a centralized batch 
mode, positions can be retrieved (up to a rigid transformation) by applying Princi¬ 
pal Component Analysis (PCA) on a so-called similarity matrix built from the rel¬ 
ative distances. In this paper, we propose a distributed on-line algorithm allowing 
each node to estimate its own position based on limited exchange of information 
in the network. Our framework encompasses the case of sporadic measurements 
and random link failures. We prove the consistency of our algorithm in the case 
of fixed sensors. Finally, we provide numerical and experimental results from both 
simulated and real data. Simulations issued to real data are conducted on a wireless 
sensor network testbed. 
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1. Introduction 


The problem of self-localization involving low-cost radio devices in WSN can 
be viewed as an example of the internet of things (loT). The evolution in the last 
50 years of the embedded systems and smart grids has contributed to enable the 
WSN integrates the emerging system of the loT. Recently, advanced applications 
to handle specific tasks require the support of networking features to design cloud- 
based architectures involving sensor nodes, computers and other remote compo¬ 
nent. Among the large range of applications, location services can be provided 
by small devices carried by persons or deployed in a given area, e.g. routing and 
querying purposes, environmental monitoring, home automation services. 

In this paper we investigate the problem of localization in wireless sensor net¬ 
works (WSN) as a particular application of principal component analysis (PCA). 
We assume that wireless sensor devices are able to obtain received signal strength 
indicator (RSSI) measurements that can be related to a ranging model depending on 
the inter-sensor distances. The multidimensional scaling mapping method (MDS- 
MAP) consists in applying PCA to a so-called similarity matrix constructed from 
the squared inter-sensor distances. Then, the sensors’ positions can be recovered 
(up to a rmid transformation) from the principal components of the similarity ma¬ 
trix yj], yD. As opposed to time difference of arrival (TDOA) and angle of arrival 
(AOA) techniques, the MDS-MAP approach allows to recover the network config¬ 
uration based on the sole RSSI, and can be used without any additional hardware 
or/and synchronization specifically devofe fo self-localization. 

MDS-MAP has been extensively studied in the literature (see Section l23] for an 
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overview). The algorithm is generally implemented in a centralized fashion. This 
requires the presence of a fusion center which gather sensors’ measurements, com¬ 
putes the similarity matrix, performs the PCA, and eventually sends the positions 
to the respective sensors. In this paper, we provide a fully distributed algorithm 
which do not require RSSI measurements to be shared. In addition, our algorithm 
can be used on-line. By on-line, we mean that the current estimates of the sen¬ 
sors’ positions are updated each time new RSSI measurements are performed, as 
opposed to batch methods which assume that measurements are collected prior to 
the localization step. Therefore, although we assume throughout the paper that the 
sensors’ positions are fixed, our algorithm has the potential to be generalized to 
moving sensors, with aim to track positions while sensors are moving. 

The paper is organized as follows. In Section|2j we provide the network and the 
observation models. We also provide a brief overview of standard self-localization 
techniques for WSN. Section [3]presents the centralized version of the MDS-MAP 
algorithm. The proposed distributed MDS-MAP algorithm is provided in Sec¬ 
tion in An additional refinement phase is also proposed in Section [5] where our 
MDS-MAP algorithm is coupled with a distributed maximum-likelihood estima¬ 
tor. In Section m numerical experiments based on both simulated and real data are 
provided. Section |7] gives some concluding remarks. 

2. The framework 

2.1. Network model 

Consider N agents {e.g. sensor nodes or other electronic devices) seeking to 
estimate their respective positions defined asj^i,-- - ,zn] where for any i, Zi G 

with p = 2 or 3. We assume that agents have only access to noisy measurements 
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of their relative RSSI values. More precisely, each agent i observes some RSSI 
measurements Pij associated with other agents j / i. Here, Pij is a random 
function of the Euclidean distance dij = \\zi — Zj\\ between nodes i and j. The 
statistical model relating RSSI values to inter-sensor distances is provided in the 
next paragraph. 

The goal is to design a distributed and on-line algorithm to enable each sensor 
node to estimate its position Zj from noisy measurements of the distances. Before 
going further in the description of the RSSI statistical model, it is worth noting that 
the localization problem is in fact ill-posed. Since the only input data are distances, 
exact positions are identifiable only up to a rigid transformation. Indeed, quanti¬ 
ties {dij)\/ij are preserved when an isometry is applied to the agents’ positions, 
i.e. rotation and translation. The problem is generally circumvented by assum¬ 
ing a minimum number of anchors or also named landmarks (sensor nodes whose 
GPS-positions are known), e.g. M = 3 or 4 when p = 2, and considering these 
prior knowledge to identify the indeterminacy. This point is further discussed in 
Section 1231 


2.2. Received signal model 

We rely on the so-called log-normal shadowing model (LNSM) to model RSSI 
measurements as a function of the inter-sensor distance OD. We define the average 
path loss PL((i) at a distance d expressed in dB as PL((i) = PLq -|- lOr/logj^Q 
where the parameters rj, do and PLq depend on the environment (see Section [^. 
Given that the distance between sensors i and j is dij, we define the RSSI between 
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i and j as a random variable Pij satisfying 


Pi,j — 


( 1 ) 


where (e^j- '■ i j) are thermal noises assumed independent with zero mean and 
variance cr^. Assume that a given agent i is provided with T independent copies 
Pjj (l),..., Pij{T) of the random variable Pij and let Pij = T~^ Ylt=i 


be the empirical average. An unbiased estimate of the squared distance dfj is given 
by 


D{i,j) 


10 


577 


( 2 ) 


cr^ In 10 

where C = . Indeed, it can be easily checked that the mean and variance 

of the unbiased estimator ^ are respectively: K[D{i,j)] = dfj and E,[{D{i,j) — 
dfj)'^] = d‘fj{C^ — 1). The construction of unbiased estimates of squared distance 
will be the basic ingredient of our distributed MDS-MAP algorithm. 


2.3. Overview of some localization techniques 

Several overview papers have been published in the last ten years dealing with 
the classification of the localization techniques (see or BSl]). In some situations, 
localization is made easier by the presence of anchor-nodes whose positions are 
assumed perfectly known. Other methods, called anchor-free, do not require the 
presence of such landmarks. 

Anchor-based methods: The classical techniques involve the resolution of a 
single unknown position of a sensor node at a time by means of RSSI values fol¬ 
lowing the LNSM coming from a fixed number of surrounding anchor nodes or 
landmarks. Since fhe sensor node only uses fhe informafion from known posi- 
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tions, its position can be expressed in absolute coordinates, i.e. anchor positions 
in GPS-coordinates. When considering a noisy scenario, several works coupled 
the classical methods (trilateration, multilateration fl] omiin-max |Q]) with a least 
squares problem. In particular, [gl] and [Tl] consider multi-hop communications 
between the sensor nodes. Other approaches focus on the statistical distribution of 
the received RSSI measurements coming from the landmarks. The goal is to con¬ 
sider a parametric model for the received signal and to apply maximum likelihood 

fl] 


estimator (MLE). Most works consider the LNSM (see for instance 
while others assume alternative statistical models (see lll2l] or ilSH l. 


or m) 


Anchor-free methods: The configuration of the network can be recovered on 
a relative coordinate system instead of the GPS absolute coordinate system. When 
distances between nodes are view as similarity metrics, the positioning problem 
is referred to multidimensional scaling (MDS). The aim is to find an embedding 
from fhe N nodes such fhaf disfances are preserved. In classical MDS iQ, Chapfer 
12] posifions are obfained by principal componenf analysis (PGA) of a x A" 
mafrix consfrucfed from fhe Euclidean disfances. If disfances are issued fo some 
noise, e.g. esfimafed from RSSI measuremenfs as Q, yD propose a MDS-MAP 
algorifhm based on fhe classical MDS problem. Indeed, fhe WSN localization 
problem is solved by enabling each sensor node fo infer all fhe esfimafed pairwise 
disfances. Alfernafive approaches wifhin fhe localizafion confexf are based on opfi- 
mizafion fechniques. In mefric MDS, posifions are obfained by fhe sfress majoriza- 


fion algorifhm SMACOE (see Uj, Chapfer 8] and jldll '). Alfernafively, semidefinife 


programming (SDP) can be used as in ilSll . 


The laffer approaches have been also addressed in a disfribufed selling wilh- 
oul fhe presence of a cenfral processing unil. A disfribufed balch version of fhe 
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SMACOF algorithm based on a round-robin communication scheme is proposed 


in il6l] . Since jig] considers the minimization of the non-convex stress function, 


the same distributed approach (batch and incremental) is presented in jlTI] but us 


ing a quadratic criterion which includes the information from the anchor nodes to 


overcome the non-convex issue. The Authors of ilSh propose a distributed imple 


mentation of their SDP-based localization algorithm. In jlSh the network is divided 


in several clusters of at least two anchor nodes and a large number of sensor nodes 
and then the SDP problem is addressed locally at each cluster. More recently, 
gossip-based algorithms have been proposed in H, H to solve the distributed 
optimization problem via Kalman filtering and gradient descent approaches. Other 
works address the distributed WSN localization problem using the multidimen¬ 


sional scaling (MDS) method based on PCA. The MDS-MAP proposed in 120 is 
later improved in |21]. In Q] each sensor node applies the MDS-MAP of iQ] 


to its local map and then the local maps are merged sequentially to recover the 

y, in 13] 


global map. Alternatively, 


and Il23ll a sparsification matrix model on the 


observations is introduced to decentralized the PCA step. 


3. Centralized MDS-MAP 

3.1. Centralized batch MDS 

Define S as fhe N x N mafrix of square relative distances i.e., S{i,j) = dfj. 
Define 'z = -^ cenfer of mass (or barycenter) of fhe agents. Upon 

noting that d‘f ^ = \\zi — zp -|- \\zj — z|p — 2(zj — z, zj — z), one has: 

S = cl^ + Ic^ - 2ZZ^ (3) 
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where 1 is the N x p matrix whose components are all equal to one, c = (||zi — 
z|p, • • • , ||z 7 v — and the ith line of matrix Z coincides with the row-vector 

Zi — z. Otherwise stated, the ith line of Z coincides with the barycentric coordi¬ 
nates of node i. Define J = 11^/A^ as the orthogonal projector onto the linear 
span of the vector 1 = (1,..., 1)^. Define = /at — J as the projector onto 
the space of vectors with zero sum, where /at is the N x N identity matrix. It is 
straightforward to verify that J±Z = Z. Thus, introducing the matrix 

M ^ -^J±SJ^ , (4) 


equation (l3]l implies that M = ZZ'^. In particular, M is symmetric, non-negative 
and has rank (at most) p. The agents’ coordinates can be recovered from M (up 
to a rigid transformation) by recovering the principal eigei^ace of M i.e., the 
vector-space spanned by the pth principal eigenvectors (see 111 Chapter 12]). 

Denote by the eigenvalues of M in decreasing order, i.e. Xi > ■ ■ ■ > 

Xj\f. In the sequel, we shall always assume that Xp > 0. Denote by 
corresponding unit-norm x 1 eigenvectors. Set Z = (\/AiUi, • • • , ^J^Up). 
Clearly M = ZZ'^ = ZZ and Z = RZ for some matrix R such that RRiF = 
I TV- Otherwise stated, Z coincides with the barycentric coordinates Z up to an 
orthogonal transformation. In particular, the fth row of matrix Z is an estimate 
of the position of the fth sensor (up to the latter transformation common to all 
sensors). In practice, matrix S is usually not perfectly known and must be replaced 
by an estimate S. This yields the Algorithm [U (see ll, Chapter 12]). 





Algorithm 1: Centralized batch MDS-MAP for localization 
Input: Noisy estimates of the square distances D{i,j) Q for all pair i,j. 

1. Compute matrix S = j))i,j=i,...,Ar. 

2. SetM = -y^SJ±. 

3. Find the eigenvectors and eigenvalues of M. 

Output: Z = (v^ui, • • • , s/^Up) 


3.2. Centralized on-line MDS 

In the previous batch Algorithm [T] measurements are made prior to the estima¬ 
tion of the coordinates. From now on, observations are not stored into the system’s 
memory: they are deleted after use. Thus, agents gather measurements of their 
relative distance with other agents and, simultaneously, estimate their position. 

3.2.1. Observation model: sparse measurements 

We introduce a collection of independent r.v.’s {Pij{n) : i, j = I, ■ ■ ■ , N, n £ 
N) such that each Pij{n) follows the LNSM described in Section 1231 At time 
n, it is possible to define an unbiased estimate Dn{i,j) the squared distance as 

-PjjM-Pro 

Dn{i,j) = ^ iri the sense that E[D„(f,y)] = d‘f-. We use the con¬ 

vention that Dn{i, i) = 0. 

Definition 1 (Sparse measurements). At each time instant n, we assume that with 
probability qij, an agent i is able to obtain an estimate of the square dis¬ 

tance with an other agent j i and makes no observation otherwise. Thus, one can 
represent the available observations as the product Sn{i,j) = Anii, j)Dn(i, j) 
where {An)n P tin i.i.d. sequence of random matrices whose components An{i,j) 
follow the Bernoulli distribution of parameter Qij. Stated otherwise, node i ob¬ 
serves the ith row of matrix An o Dn at time n where o stands for the Hadamard 
product. 
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Lemma 1. Assume qij > Ofor all pairs i,j. Set W := [qij^]^j=i and let An, Sn 
be defined as above. The matrix 


Sn = WoAnODn 


(5) 


is an unbiased estimate of S i.e., IE[S^n] = S. 

Proof. Each entry of matrix Sn, Snii,j), is equal to l/qijAn{i,j)Dn{i,j). As 
the random variables An{i,j) and Dn{i,j) are independent, by the above defini¬ 
tion of Dn and E[A„(z, j)] = qij, then E[S'n(i, j)] = df -. □ 

As a consequence of Lemma [T] an unbiased estimate of M defined in (|4ll is 
simply obtained by Mn = —^J±SnJ±. 


3.2.2. Oja’s algorithm for the localization problem 

When dealing with random matrices Mn having a given expectation M, the 


principal eigenspace of M can be recovered by the Oja’s algorithm 12411 . The latter 


consists in recursively defining a sequence Un of N x p matrices, which stand for 
the estimate at time n of the p principal unit-eigenvectors of M. The iterations as 


firstly introduced in 12411 are given by: 


Un — U n—1 ~\~ 'Yn {^MlnU n—1 Un(U n—1-^ nU , 


( 6 ) 


where > 0 is a step size. Note t 


rat in practice, the algorithm is likely to suffer 


from numerical instabilities. In i25h . a renormalization step is introduced to avoid 


unstabilities. As this approach seems difficult to generalize in a distributed context. 
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it is more adequate in our context to introduce a reprojection step in ® of the form 


Un — IIjc \Un-1 +7n nUn-1 “ Un{U'^_iMnUn-l))] , 


where IIjc is a projector onto an arbitrarily large convex compact set % chosen 
large enough to include all matrices whose columns have unit-norm. Typically, we 
set X = [—a, aY x • • • x [—a, aY where a > 1. 

In order to obtain an estimate of the sensors positions, we also need to estimate 
the principal eigenvalues in addition to the eigenvectors. Let ^ denote the kth 
column of matrix Un- Define the quantity Xn,k recursively by: 


^n,k ^n—l,k T In {X^n—l,k-^^n'^n—l,k l,fc) 


(V) 


The convergence properties of Oja’s algorithm are studied in details in 1241] and 


I 25h . Finally, according to step 3 of the batch Algorithm [T] the estimated barycen- 
tric coordinates are obtained as: 


Zn — 


• • • ) \/ Xn^pUn^p'j 


( 8 ) 


The combination of Equations ® O and ([8]l provides an on-line for MDS-MAP 
algorithm. However, the computation of matrix at each step as well as the 
matrix products in ® require a full amount of centralization. 
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4. Distributed on-line MDS-MAP 


4.1. Communication model 


It is clear from the previous section that an unbiased estimate of matrix M 
is the first step needed to estimate the sought eigenspace. In the centralized set¬ 
ting, this estimate was given by matrix = —^J±SnJ±- As made clear by 
the observation model (in Definition [B, each node i observes the ith row of ma¬ 
trix Sn- As a consequence, node i has access to the ith row-average Sn{i) — 
Sn{i,j)- This means that matrix SnJ± can be obtained with no need to 
further exchange of information in the network. On the other hand, J±SnJ± re¬ 
quires to compute the per-column averages of matrix SnJ±, i-e. Yj 
for all i. This task is difficult in a distributed setting, as it would require that all 
nodes share all their observations at any time. A similar obstacle happens in Oja’s 
algorithm when computing matrix products, e.g. MnUn-i in (0). To circumvent 
the above difficulties, we introduce the following sparse asynchr'onous communi¬ 
cation framework. In order to derive an unbiased estimate of M, let us first remark 
that for all i, j, 




d‘^{i) + _ dlj + S 

2 2 


(9) 


where we set d'^{i) = ^ Yk ^ ~ Note that the terms d'fj and 

d'^{i) can be estimated by Snii,j) and Sn{i) respectively. However, additional 
communication is needed to estimate <5 since it corresponds to the average value 
over all square distances. We define 




Sn{i) + Sn{j) SniiJ) + Sn{i) 


( 10 ) 
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where Sn{i) is a quantity that we will define in the sequel, and which represents 
the estimate of <5 at the agent n. 

We are now faced with two problems. First, we must construct (5„(i) as an 
unbiased estimate of 6. Second, we need to avoid the computation of Mn{i,j) for 
all pairs i,j, but only to some of them. In order to provide an answer to these prob¬ 
lems, we introduce the notion of asynchronous transmission sequence. Formally, 


Definition 2 (Asynchronous Transmission Sequence). Let qbe a real number such 
that 0 < q < 1. We say that the sequence of random vectors = {tn,Qn,i ■ 
i G {1, • ■ ■ ) N},n G N) is an Asynchronous Transmission Sequence (ATS) if: i) 
ail variables {tn,Qn,i)i,n tire independent, ii) tn is uniformly distributed on the 
set {1, • • • , A^}, iii) Vi in, Qn,i A a Bernoulli variable with parameter q i.e., 
F[<5n,i = 1] = g and iv) Qn,,^ = 0. 


Let {Tn)n denote an ATS defined as above. At time n, we assume that a given 
node in G {1,..., A^} wakes up and transmits its local row-average Sn{in) to 
other nodes. All nodes i such that Qn,i = 1 are supposed to receive the message. 
For any i, we set: 


^n(i) 


Sn{i) 

N 


+ 


Sn fn)Qr 


( 11 ) 


The following Lemma is a consequence of Definition |2] along with Lemma [U 
and equation dUl. 


Lemma 2. Assume that {Tn)n A an ATS independent of(Sn)n- Let (Mn)n be the 
sequence of matrices defined by AlOh . Then, E[AT„] = M. 


Proof By Lemma [T] the expectation of terms Sn{i), Sn{j) and Sn{i,j) are re¬ 
spectively d?‘{i), d?{j) and Moreover, by Definition [2] the expectation of the 
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random term Sn{i) is equal to 


^ j¥=i t=l 

which coincides with <5. Then, the expectation of each entry of the matrix iW„ 
in (fTOl) is equal to the corresponding M (i, j) defined in (|9ll. □ 

4.2. Preliminaries: constructing unbiased estimates 

As we now obtain a distributed and unbiased estimate of M, the remaining 
task is to adapt accordingly the Oja’s algorithm ([61). In this paragraph, we provide 
the main ideas behind the construction of our algorithm. 

Assume that we are given a current estimate Un-i at time n, under the form of 
a xp matrix. Assume also that for each i, the fth row of Un-i is a variable which 
is physically handled by node i. We denote by the fth row of Un-i- 

Looking at ® in more details, Oja’s algorithm requires the evaluation of inter¬ 
mediate values, as unbiased estimates of MUn-i and n-i- 

We consider the previous ATS {Tn)n involved in (fTOl ). We assume that the 
active node in (i.e., the one which transmits Sn{tn)) is also able to transmit its 
local estimate Un-i{tn) at same time. Thus, with probability jj, node in sends its 
former estimate Un-i{tn) and Sn{tn) to all nodes i such that Qn^i = 1. Then, all 
nodes compute: 

— iV — 

Yn{i) = Mn{i,i)Un-l{i) H- Un-l{tn)Mn{i,in)Qn,i (12) 

q 

As it will be made clear below, the N xp matrix Yn whose zth row coincides with 
Yn{i) can be interpreted as an unbiased estimate of MUn-i- 
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Now we introduce the distributed version of the second term nUn-i- 

Consider a second ATS {T!^)n independent of {Tn)n- At time n, node wakes up 
uniformly random and broadcasts the product to other nodes. 

Receiving nodes are those i’s for which Q'^ i = 1- Then, all nodes are able to 
compute the estimate p x p matrix as follows: 

= Un-l{i) Yn{i) —Un-l{Ln) nil'n)Qn,i ■ (1^) 

Lemma 3. Let {Tn)n and {T!^)n be two independent ATS. For any n, denote by 
the a-field generated by {Tk)k<n, iU) k<nj (-^A:)fc<n Let (Un)n 

a 3'n-kneasurable N xp random matrix and let Y^, be defined as above. Then, 

= MUn-i and E[A„(i)|T„_i] = Ul_^MUn-i ■ 

Under Lemma\I}^and Definition^ the random sequences Yn{i) and A„(i) are 
unbiased estimates M{i, j)Un-i{j) and U'^_iMU n-i respectively given 

Un-l. 

Proof. For each i, we obtain 

E[y„(i)|J„_l] = M{ifi)Un-l{i) + -j^Y.M{i,j)Un-l{j) 

3 
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and 


]E[A„(i)|J„_i] = Un-i{i)'^nYn{mn-i] + -^ J7„-i(j)^E[y„(j)13-„_i]g 

= E E ^n-l(i)^M(i, j)Un-l{j) 
i j 

which corresponds with the square matrix Ujl^_iMUn-i- □ 

4.2.1. Main algorithm 

We are now ready to state the main algorithm. The algorithm generates itera¬ 
tively and for any node i two variables Un{i) and A„(i), according to: 

Unii) = Un-lii) + In C^n(i) “ Un-l(i}An(i)} (14) 

An(i) — Aj^_i(z) + yn,(diag(A,2(^)) Afi_i(i)). (15) 

For the same reasons as before, it is important in practice to introduce a projection 
step rtuc in (fT4l) to avoid numerical unstabilities. Finally, as in ([8]l, each sensor i 
obtains its estimate position Zn(i) by: 

Zn(i) — Aj^ 1 , ■ ■ ■ ) Y^Anip(i)rifip(i)^ (16) 

where we set Unii) = iun,i{i), ■ ■ ■ ,Un,p{i)). 

The proposed algorithm (fT4]) -(fT^ is summarized in Algorithmic below. Note 
that, at each iteration time n, only two communications are performed by two 
randomly selected nodes issued to the ATS’s T„ and T^. 
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Algorithm 2: Distributed on-line MDS-MAP for localization (doMDS) 
Update: At each time n = 1, 2,... 

[Measures]: each sensor node i, do: 

Makes sparse measurements of their RSSI to obtain {Dn{i,j))j for some j 
such that An{i,j) = 1 (Definition [T]). Set 


Sn{iJ) 


Qi/Dnii,j) ifAn{i,j) = l 

0 otherwise 


and set 5^1 (z) jy >S'j^(z, j). 

[Communication step]: 

A randomly selected node z„ wakes up, then 

i) The node Ln randomly selected broadcasts Un-i{i-n) and Sn{tn) to 
nodes i such that Qn^i = 1 - 

ii) Each node i computes Yn{i) by ([12]). 

Hi) A node z'„ randomly selected broadcasts Unii^n) to 
nodes i such that Q'^ j = 1 - 

iv) Each node i updates Un{i) by (IT3]) - ([T4]) and Zn{i) by ([T^ . 


4.3. Convergence analysis 

We make the following assumptions. The sequence ( 7 n)n is positive and satis¬ 
fies 

^ 7 n = +00 and ^ 7 ^ < oo . 

n n 

Moreover we make the assumption that the sequence Un remains a.s. in a fixed 
compact set %. It must be emphasized that this assumption is difficult to check in 
practice. As mentioned above, stability can be enforced by means of a projection 
step onto %. 

Proposition 1. For any U G set h{U) = MU - UU'^MU. Let Un be 

defined by Aldt . There exists a random sequence such that, almost surely (a.s.), 
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Yin Inin Converges and 


Un — U n—l + 'ynh{Un—l) Inin ■ ( 17 ) 


The proof is provided in the Appendix. We are now in position to state the 
main convergence result. 

Theorem 4. For any k = 1, • • • ,p, the kth column Un,k of Un converges to an 
eigenvector of M with unit-norm. Moreover, for each node i, Xn,k{'i) converges to 
the corresponding eigenvalue. 


The proof is provided in the Appendix. 

Note that Theorem |4] might seem incomplete in some respect: one indeed ex¬ 
pects that the sequence Un converges to the principal eigenspace of M. Instead, 
Theorem m only guarantees that one recovers some eigenspace of M. Undesired 
limit points can be theoretically avoided by introducing an arbitrary small Gaus¬ 
sian noise inside the parenthesis of the left hand side of (fTdl) (see Chapter 4 in 


ll26n i. These so-called avoidance of traps techniques are however out of the scope 


of this paper, and numerical results indicate that the principal eigenspace is indeed 
recovered in practical situations. 


5. Position refinement: distributed maximum likelihood estimator 


In the context of WSN localization, a refinement phase is in general added 
(see |9], iQl, 121] or 116]). It is usually based on the statistical model relating the 
observed RSSI values to the unknown positions, the latter being estimated in the 
maximum likelihood sense. The objective is twofolds. First, maximum likelihood 
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estimation improves the estimation accuracy obtained by the MDS-MAP approach. 
Second, as the MDS-MAP only identifies positions up to a rigid transformation, it 
allows to eliminate the residual ambiguity by using anchor nodes, provided that 
such anchors exist. 

In this section, we provide a distributed algorithm in order to locally maxi¬ 
mize the likelihood. It is worth noting that the likelihood function is generally 
non-convex. Thus, one cannot expect that a standard gradient ascent provides the 
maximum likelihood estimator regardless from the initialization. For this reason, a 
preliminary phase such as the proposed doMDS algorithm is essential as an initial 
coarse estimate, and the algorithm depicted below is used merely as a fine search 
in the neighborhood of the doMDS output. 

5.1. Likelihood function 

Consider a connected graph G = {V,E) where V = {l,...,A^}is the set of 
agents and £" is a set of non-directed edges. In this paragraph, we allow for the pres¬ 
ence of anchor nodes. We let A C {1,... , A^} be the set of anchor nodes i.e. for 
each k £ A, the position of node k is assumed to be known. Unknown parame¬ 
ters thus reduce to set of coordinates z = (zj : z S A) where A = U\A. We denote 
by jVi the neighbors of i which belong to A and by .Jfi the neighbors of i which are 
anchors. We note z_ 4 < = (zj : j G Ai U {z}). For a connected pair of nodes {z, j}, 
we let Pij{n) in G N) be an i.i.d. sequence following the LNSM model of Sec¬ 
tion |2]2] Equivalently, the quantity iij{n) = follows a normal dis¬ 

tribution with mean log^g di,j and variance since h,j{n) = log^o dij + ^ 
by using ([Til. Based on the observations {ii,j{n) : z ~ j) at a given time n, the 
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likelihood associated with the unknown sensors’ positions can be decomposed as 


N 

i=l 

where 


E 


fi{z._A^i)= > fij(n)-logio 


\Zi - z. 


+ X] - logio ll^:* - Zfc 


5.2. The algorithm: on-line gossip-based implementation 


Following the idea of 12711 (see also 12811 and reference therein), we propose a 


distributed consensus-based implementation consisting on local computations and 
random communications among the sensor nodes. The algorithm is given below. 
The convergence proof is omitted due to the lack of space but follows from the 


same arguments as ll28h . 


Algorithm 3: Distributed on-line MLE (doMLE) 

Update: at each time n = 1, 2,... 

[Local step] each node i obtains and {Pik{n)}\tk&jii- 

Each sensor i computes a temporary estimate of its position’s set: 


[Gossip step] two uniformly random selected nodes i ~ y in A exchange 
their temporary estimated positions and average their values according to: 




^^- 


Otherwise, Ml ^ jVj or m i,j, set = z^.^^n{l)- 


Algorithm [3] uses a standard pairwise averaging between nodes. We note that 
more involved gossip protocols have been proposed, we mention for instance broad- 
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cast and push-sum protocols (see 12911 and OOl] '). Although theoretically possible, 
such an extension of Algorithm |3]is however beyond the scope of this paper. 


6. Numerical results 


We consider the same network configuration corresponding on the set of = 
50 sensor nodes selected from the FIT loT-LAB 1^ platform at Rennes. Sensor 
nodes are located within a 5 x 9m^ area, i.e. p = 2. Six sensors of the 50 were 
set as anchor nodes (or landmarks). We compare the performance of our proposed 
distributed on-line MDS (doMDS) to other existing algorithms. We consider the 


distributed batch MDS lllal (dwMDS) and the classical centralized methods such 
as: multilateration 16] (MC), min-max |7], Algorithm [T]in Section iTT] (batch MDS) 
and the Oja’s algoritm ©-dVll described in Section lT2l The three iterative algo¬ 
rithms (Oja’s, dwMDS and doMDS) are initialized by randomly chosen positions 
in 5 X 9 m^. 


6.1. Simulated data 

First, we show the results from simulated data drawn according to the obser¬ 
vation model defined in Section 13.21 In order to compare our proposed algorithm 


with the distributed MDS proposed by jlhll . we set the same environmental context 


in which a/rj = 1.7. Figure [T] displays the comparison of the root-mean square 
error (RMSE) when running Algorithm |2] over 300 independent runs of the es¬ 
timated positions when considering different communication parameters: {qij)ij 
(the Bernoullis related to the observation model dS])) and q (the Bernoullis related 


^FIT loT-LAB https : / / www .iot-lab.info/ 
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to the ATS in Definition O. Since the variance of the error sequence is upper 
bounded by the minimum probability value in (IA.3I )- (IA.5I ). we observe from Fig¬ 
ure [U a trade-off between the accuracy and the number of communications as the 
RMSE decreases faster when the probability q is closer to 1. 

Figure [2a] shows the comparison of the localization RMSE over 300 indepen¬ 
dent runs of the overall estimated positions when considering the three iterative 
methods: the centralized Oja’s ©-([T]), the dwMDS of jig] and our proposed Al¬ 
gorithm |2] The estimated positions after 1000 iterations of the three iterative al¬ 
gorithms are reported in Eigure |2| Note that, the result in Eigure [2^ requires at 
least twice the number of communications compared to the results both on-line 
Oja’s approaches. Positions close to the barycentric of the network tend to be more 
accurate than positions on the surrounding area for the three cases. Nevertheless, 


Eigures [2b| and |2d| show these outer positions better preserved than ilol . Indeed, 


our distributed and asynchronous Oja’s algorithm achieves in general better accu¬ 
racy (around the 65% of positions) except for the third part of nodes which are 
located around the network’s boundary, e.g. nodes 11 or 36 — 37 for instance (see 
squared nodes in Eigure l2dl). 


6.2. Real data: FIT loT-LAB platform of wireless sensor nodes 
6.2.1. Platform description 

In order to obtain real RSSI values we make use of the EIT loT-EAB platform 
deployed at Rennes (Erance). The 256 WSN430 open nodeavailable at the plat¬ 
form are issued to the standard ZigBee IEEE 802.15.4 operating at 2.4 GHz. The 

^Seethe technical specifications of WSN430 sensors https :/ / git hub. com/iot- lab/iot-lab/wiki/Hardware_Wsn4 
and CC2420 transceivers involve in our campaigns: http; / /www. ti . com/lit/ds/symlink;/cc2 42 0 . pdf 
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sensor nodes are located in two storage rooms of size 6x15 containing differ¬ 
ent objects. They are placed at the ceil which is 1.9 m height from the floor in a grid 
organization. Through of our user profile created in the FIT loT-LAB’s website, 
we run remotely several experiments involving the 50 selected sensor nodes within 
5 X 9m^. All real data used in this section can be found in 3- The environment 
parameters issued to the LNSM ([Til are: = 28.16 dB, PLq = —61.71 dB and 

r/ = 2.44. We set qij = 0.8 Vi, j, q = 0.85 and 7 „ = for Algorithm[2l 

6.2.2. Performance comparison 

We compare the same algorithms considered in Section IhTI by setting the esti¬ 
mated positions obtained from each algorithm to the initialization of Algorithm [3l 
Table[T]shows the RMSE values before and after the refinement phase. In addition, 
we include the ratio of the accuracy improvement considering the RMSE values af¬ 
ter and before applying the distributed MEE and the ratio regarding the number of 
positions over the total N that are improved. The best performances are achieved 
by min-max, dwMDS and doMDS in terms of minimum RMSE value over the N 
estimated positions. Nevertheless, the highest improvement is obtained with the 
proposed doMDS since the RMSE before the refinement phase was higher than 
the values from min-max and dwMDS which do not experiment a considerable de¬ 
crease. In general, the refinement Algorithm |2] improves almost all the positions 
for each method and especially the anchor-free methods based on the MDS ap¬ 
proach. Indeed, the highest values are those from the distributed versions which 
may exploit in advantage the local knowledge of each sensor node. 

"^Database available at G. Mortal personal website http; / /per so . telecom-paristech. f r / -morralad/ 
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7. Conclusion 


This paper introduced a novel algorithm based on Oja’s algorithm for self¬ 
localization in wireless sensor networks. Our algorithm is based on a distributed 
PCA of a similarity matrix which is learned on-line. Almost sure convergence of 
the method is demonstrated in the context of vanishing step size. The algorithm can 
be coupled with a distributed maximum likelihood estimator to refine the sensors 
positions if needed. Numerical results have been conducted on both simulated and 
real data on a WSN testbed. Although we focused on fixed sensors posifions, fhe 
on-line nafure of fhe algorifhm makes if suifable for use in dynamic environmenfs 
where one seek fo frack fhe posifion of moving sensors. 

Appendix A. Proof of Proposition [T] 

Sef for each i, M{i, j)U and 

i^{i) = {Yn{i) - {MUn-l)i) + (A.l) 

Then, fhe sequence generafed by each sensor node i is wriffen as: 

Un{i) = Un-l{i)+ln {{MUn-l\ - Un-l{i){Ul_^MUn-l)) + 7nCn(i) 

Denofe by the cr-algebra generated by all random variables defined up fo fime n. 
Using Lemmas 1 1 121 and [3l if is immediafe fo check fhaf = 0 and fhus 
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the sequence Ylk<n Ik^k is 3'n-adapted martingale. We estimate 


E[||en(i)f |5'n-i] < E[||y„(i)||2|J„_i] + ||J7„_i(i)||2E[||A„(z)||2|5-„_i] 

+ 2||J7„_i(z)||E[||y„(i)A„(i)|||?-„_i]. (A.2) 


The first term on the right hand side (RHS) of (IA.2I) can be expanded as: 

E[||l^„(z)f |T„_l] < n\Mn{hi)\^]\\Un-l{i)f 

+ -Y,n\Mn{i,j)n\Un-l{i)f 

+ 2Y,¥.[Mn{hi)Mn{hj)]\\Un-i{i)f. (A.3) 

j¥=i 

Upon noting that for any i,j K[Sn{i,j)‘^] = and Un-i lies in a fixed 

compacf sef, fhere exisfs a consfanf K' such fhaf E[||l^„(f)|p|3'n_i] < K' for all 
n depending on N, qmin = niiiij j qij, C defined in Q and maxj j dfj such fhaf 
E[|Af„(i, j)p] < K for some consfanf K. The second ferm on fhe RHS of (IA.2I ) 
can be handled similarly: 

E[||A„(f)f|T„_i] < E[||r4z)||2|T„_i]||[/„_i(f)f + (^^E[||r4j)||2|T„_i] 
+ 2^E[l^„(i)l^„(j)|T„_i])||f7„_i(i)f < A" 

(A.4) 


25 






for some constant K”. Finally, 


+ ^E[r4z)y„(j)|J„_i]||J„_i(j)|| (A.5) 

j¥=i 

is uniformly bounded as well. Therefore, we have shown that a.s. 


supE[||^„(i)|p|g“„_i] < oo 

n 

Since < oo, it follows that X^„7n]E[||^n(i)|P|3'n-i] < oo a.s. By Doob’s 

Theorem, the martingale IkCki^) converges almost surely to some random 

variable finite almost everywhere. This completes the proof. 


Appendix B. Proof of Theorem |4] 

Consider the following Lyapunov function V : \ {0} —^ R+: 


V{U) 


e\\u\\^ 

U'^MU ’ 


(B.l) 


The following properties hold: 

i) lim||; 7 ||^oo ’L([/) = +00 and the gradient is W{U) = —2-^^^^h{U). 

ii) {V{U), h{U)) < 0 and the equality holds iff {[/ e \ h{U) = 0}. 


The proof is an immediate consequence of Proposition [H the existence of (IB.Ill 


along with Theorem 2 of Olll . Sequence Un converges a.s. to the roots of h. The 


latter roots are characterized in 12411 . In particular, h{U) = 0 implies that each 


column of U is an unit-norm eigenvector of M. 
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Figure 1: RMSE as a function of nN from the two estimated eigenvectors Un,i and 
Un ,2 when considering the noiseless and noisy case and for different values of q. 


Method 

MC 

min-max 

MBS 

Oja 

dwMDS 

doMDS 

Before refinement 

1.87 

0.8 

1.98 

2.18 

0.86 

1.56 

After refinement 

1.05 

0.54 

1.39 

1.37 

0.6 

0.51 

Improvement (%) 

44 

32 

30 

28 

30 

78 

Positions improved (%) 

75 

71 

80 

80 

82 

86 


Table 1: RMSE averaged over the 44 estimated positions considering real data. 


31 



















(a) RMSE as a function of nN from the estimated positions (Z„(l),..., Zn{N)). 
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(b) Oja’s algorithm (|6l)-(|7]i. (c) dwMDS ifl^ . 


(d) Algorithm|2](doMDS). 


Figure 2: Estimated positions after 1000 iterations. Markers (♦) correspond to the 
estimated values while markers (O) to the true positions. Squared positions (□) in 
d) highlight worse accuracy compared to b). 
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